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Introduction 


GAMNAS (Geometric and Material Nonlinear Analysis of Structures) is a 
two-dimensional finite element stress analysis program. The program was 
developed to support fracture mechanics studies of debonding and delamination 
(refs. 1-3). Options include linear, geometric nonlinear, material nonlinear, 
and combined geometric and material nonlinear analysis. 

The purpose of this manual is to document the theoretical basis of GAMNAS 
and to provide instruction in the use of the program. Details of the program 
organization and logic are presented in order to guide the user who needs to 
modify the code to meet some special need. Familiarity with linear finite 
element analysis is assumed. 

First, theoretical aspects of GAMNAS are presented. Then program speci- 
fications, such as allowable problem size, are given. Next the program orga- 
nization is described. Finally, the required input data is described. Brief 
descriptions of the subroutines and major program variables are given in 
Appendix A. Appendix B gives input data and results for several sample 
problems. Appendix C briefly discusses error messages and possible debug 
strategies. 

Successful use of any finite element program depends largely on the 
ability of the analyst to qualitatively predict the response of a configura- 
tion before attempting detailed finite element analysis. This insight will 
generally be based on experience and possibly some strength of materials 
arguments. Also, very coarse finite element models can be useful. Insight is 
particularly important for nonlinear analysis, in which questions of conver- 
gence, uniqueness of solution, and solution strategy must be addressed. 

Hence, the user should become thoroughly familiar with the theoretical basis 
of GAMNAS and then gain experience by analyzing a variety of simple configura- 
tions before attempting to analyze a complex configuration. 
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Nomenclature 


[B] 

[D*] 

E 

F 

G 

G I’ G II 

I 

[K] 

[KoMRp] 

M 

P-,P- 
x y 

[R] 

[R] 

[T] 


u,v 


u,v 


V 


x,y 

x,y 

3 

Aa 
{ A 6} 
{ Ae} 
{Ai} 
{Ae} 
{ 6 } 
{ 6 } 


incremental strain-displacement matrix 

elasto-plastic constitutive matrix 

Young's modulus 

yield surface function 

total strain-energy-release rate 

mode I and mode II components of strain-energy-release rate 
moment of inertia 

transformed global stiffness matrix 

linear and tangential stiffness matrices, respectively 
moment 

forces transmitted through crack tip in the x and y directions 
applied load vector 
transformed applied load vector 
transformation matrix 

displacements in x and y directions, respectively 
displacements in x and y directions, respectively 
volume 

rectangular Cartesian coordinates 
rotated rectangular Cartesian coordinates 

fraction of strain increment required to reach yield surface 
virtual crack closure length 
increment to nodal displacement vector 
strain increment 

strain increment required to reach yield surface 
strain increment after reaching yield surface 
nodal displacement vector 
transformed nodal displacement vector 
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£ , £ , £ 
x y xy 


{€} 

U, r 

X 

0 c 
er 


ys 

{a} 


{a o> 

{*> 


normal strains in x and y directions and shear strain in xy 

plane, respectively 

strain vector 

plastic strain increment 

proportionality constant 

2 2 2 

effective stress, equal to (o + a + a “ oo - a a - 

v x y z xy yz 

. O 2 , 1/2 

a a + 3a 1 
x z xy J 

uniaxial yield stress 
stress vector 

stress vector before strain increment 
stress vector after strain increment 
residual force vector 


Theory 

Governing Equations 

This section outlines the theoretical basis for the GAMNAS computer 
code. First, geometric and material nonlinearity are discussed in general. 
Then application of the displacement based finite element method to nonlinear 
problems are discussed. The description given here follows closely that given 
in refs. 4 and 5, where details may be found. 

Herein, geometric nonlinear analysis refers to an analysis which calcu- 
lates strains using the Lagrangian nonlinear strain-displacement relations, 
eqns. (1) 


_ 3u . 1 | 
e x 8x 2 1 

/9u\ 2 /9v\ 2 1 
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The second-order terms in eqns. (1) account for finite rotations* However > 
the strains are still assumed to be infinitesimal. 

For material nonlinear analysis, the nonlinear relationship between 
stress and strain is defined incrementally , eqn. (2) 


d{a} = [D*] d{ £ } (2) 

The nonlinear elasto-plastic constitutive matrix [D*] is a function of the 
assumed yield surface and flow rule and the current stress state. GAMNAS uses 
the von Mises yield surface, eqn. (3) and a flow rule based on the normality 
principle, eqn. (4). 



a o - a a 
x y y z 


a a 
x z 


30 Z ) 

xy ' 


1/2 


- 0 


ys 


(3) 


d{£} p 



(4) 


The instantaneous uniaxial yield stress, a , in eqn. (3) is a function of the 

ys 

strain history and the specified uniaxial stress-strain curve. Three types of 
uniaxial stress-strain curves can be specified: elasto-plastic, bilinear, and 

Ramberg-Osgood. These are shown schematically in fig. 1. GAMNAS can only 
analyze plastic deformation of isotropic materials. 

Application of the finite element method to nonlinear problems is very 
similar to that for linear problems. In both cases a system of equations is 
derived which expresses the equilibrium of internally generated forces in a 
body with externally applied forces, eqn. (5) 




f [B] T {a} dV - {R} = 0 
VOL 


(5) 


In eqn. (5) {if;} , {a}, and {R} are the residual force, stress, and applied load 
vectors, respectively. The integral is the vector of internally generated 


f* 
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forces. The matrix [B] is the incremental strain-displacement matrix, as 
defined by eqn. (6) 

d{e} = [B] d { 6} (6) 

where [6] is the nodal displacement vector, i.e. , a list of u and v 
displacements at the nodes. For linear problems eqn. (5) is a linear set of 
equations with unknowns (6). 

For geometrically nonlinear problems eqns. (1) are used with eqn. (6) to 
derive [B ] • The matrix [B] is found to vary linearly with {6}, as is expected 
from the quadratic form of eqn. (1). The stresses {a} are linearly related to 
the strains, which vary quadratically with (6). Hence, eqn. (5) is a set of 
cubic equations in {6}. 

For elas to-plastic problems, the matrix [B] is independent of {6}, but 
the relationship between {6} and {a} is a complicated nonlinear function. 
Furthermore, the relationship between {6} and {a} is path (i.e., history) 
dependent. Hence, the solution of eqn. (1) for a desired load level is 
obtained by dividing the total load into a series of small load increments. 

For each load increment, the relationship between stress and strain is deter- 
mined from eqn. (2). 

For combined geometric and material nonlinearity, the nonlinear relation- 
ships for each are simply used together. 

Iterative Solution 

The governing equations, eqn. (5), are solved iteratively using modified 
Newton-Raphson methods (ref. 4). The basic Newton-Raphson method for the 
first load step is outlined below. 

1. Obtain a linear solution using the linear stiffness matrix K q : 

<«„> ■ 1K o fl < R > 
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2. Calculate residuals {^} with eqn. (5) 

3. Check for convergence. Stop if {^} is sufficiently small. 

4. Calculate tangential stiffness matrix, [K T ] 

(The tangential stiffness matrix is defined by the equation 
[K^,] { A6} - {A*}.) 

5. Solve for correction to displacements 
{AS} = - [K x _1 ] {*} 

6. Update displacements: 6 = 6 + A6 

7. Go to step 2. 

If multiple load steps are used, only step 1 changes. After obtaining a 
converged solution for load step "i", the linear solution (i.e., the new 
step (1)) for the next load step is 

(S}. +1 = {6}. + [K t ]~ 1 (AR}. +1 (7) 

where {AR} . , , is the load increment, 
l+l 

Different versions of the Newton-Raphson technique described above were 
used for geometric nonlinear, material nonlinear, and combined nonlinear 
analysis in GAMNAS. The main differences are in the way the tangential stiff- 
ness matrix, [Krj, ] , is approximated. For geometric nonlinear analysis [K-p] is 
updated every "NCYCLE" iterations, where NCYCLE is an input parameter. For 
material nonlinear analysis [K T ] is approximated by the linear stiffness 
matrix [K Q ] for all iterations. For combined geometric and material nonlinear 
analysis [ICp] is updated every "NCYCLE" iterations, but the linear stress- 
strain relations are used in calculating [K^]. For combined nonlinear 
analysis the solution for each load increment begins with obtaining a con- 
verged solution in which no additional yielding is allowed. After obtaining 
this "transition" solution, iterations begin in which both geometric and 
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material nonlinear effects are included. This procedure reduces spurious 
material yielding which can be an artifact of iterative solution procedures. 
This procedure will be discussed further in the discussion of the flowchart 
for the subroutine ITERATE. 

Strain Energy Release Rates 

GAMNAS has the option to calculate Mode I and Mode II strain energy 
release rates. Strain energy release rates are calculated using a virtual 
crack extension technique similar to that reported in ref. 6. This technique 
uses the forces transmitted across the crack tip and the relative displace- 
ments just ahead of the crack tip to determine the energy release rate. For 
geometrically nonlinear problems the forces and displacements are transformed 
to the local rotated coordinate system, as shown in fig. 2. Figure 2 also 
shows the equations used to calculate Gj and Gj^. The strain energy 
release rate calculation is valid for linear and geometrically nonlinear 
analysis only . The program assumes the mesh around the crack tip is rectangu- 
lar and that the crack is initially parallel to the x-axis. Near the crack 
tip the mesh must be symmetrical about the crack tip . 


Boundary Conditions 

The following boundary conditions can be prescribed in GAMNAS: 

1 . Nodal loads 

2. Specified displacements 

3. Equivalence of two or more displacements, e.g. , 6^ = S_. 

4. Equivalence of one displacement and the negative of another 


displacement, e.g., 6 ^ = -6. 


^ |* 

To prescribe a displacement 6^ = the diagonal terra of the i c equation is 
replaced by a large number, 10^®, and the "load" terra for the i equation is 


set to 10 


30 


6 . To impose a multi-point constraint, i.e., 
o 
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^ or 6 ^ = - 6 j , the displacement and load vectors and the stiffness 

matrix are transformed (ref. 10). The transformation is best explained by 
example. Consider the linear system [K] { 6 } = {R}. Assume there are four 
nodal displacements. To impose the condition 6 ^ = 63 a new displacement 
vector { 6 } is defined such that 


{ 6 } = [T] { 6 } 



The new stiffness matrix [K] and load vector {R} are 
[K] = [T] T [K] [T] 

{R> = [T] T {R} (9) 

The new governing equations are [K] {5} = {R}. Note that ^ - S 3 . 

Hence, to impose the condition 6 ^ = 63 , we need simply impose the condition 



When multi-point constraints are specified, the bandwidth generally 
increases. The increase in bandwidth depends on the node numbering scheme. 
Hence, the multi-point constraints should be considered when selecting the 
node numbering scheme. 

Elements 

GAMNAS uses the four-node isoparametric quadrilateral. This element is 
well known to perform poorly in modeling bending type deformation when exact 
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integration is used. But the performance can be dramatically improved by 
using selective reduced integration. References 7 and 8 describe the proce- 
dure for linear problems. Reference 9 describes the procedure for geometri- 
cally nonlinear problems. The user can specify either full or selective 
reduced integration in the program. 

Program Specifications 

GAMNAS is written in Prime’s extended version of FORTRAN 77. Gore 
requirements are 604,000 16 bit words and compilation time is approximately 2 
minutes on the Prime 750. Execution times will vary greatly depending on the 
particular finite-element model. The current maximum allowable values of the 
major parameters are given in the description of the input data. An in-core 
equation solver is used. Hence, the maximum problem size is limited by the 
memory of the computer being used. 

Most of the core requirements are for holding the global stiffness 
matrix, "SN." The matrix SN is dimensioned (1300, 70), which permits 1300 
degrees of freedom (650 nodes) and a bandwidth of 70. GAMNAS can be quickly 
modified using a text editor to change the maximum bandwidth and number of 
nodes. The required changes and the order the changes should be made are 
listed below: 

1) Change the string ”(1300, 70 M to "(XXX, YYY" everywhere, where XXX and 
YYY are the new number of rows and columns, respectively. 

2) Change the string "(1300" to "(XXX" everywhere, where XXX is the new 
number of rows. 

3) In subroutine INITIAL change the following two lines: 

MRANK = 1300 + change 1300 to XXX 

MIBW = 70 «- change to 70 to YYY 
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where XXX and YYY are the new number of rows and columns in SN, 
respective ly. 


Program Organization 

In this section the flow of GAMNAS is described. Flowcharts are given 
for the more complicated routines: the main program, ITERATE, and STRSCAL. 
Very brief description of the subroutines and the major program variables are 
given in Appendix A. 

An annotated flowchart for the main program is shown in Figure 3. Only 
one proportional load vector is input. The different load numbers (LOADNUM) 
refer to the scale factor by which the load vector is multiplied. For each 
new load, a linear incremental solution is obtained in the main program before 
calling ITERATE to obtain the nonlinear incremental solutions. the linear 
solution for the first load step and all nonlinear solutions are output. 

Figure 4 shows a flowchart for the subroutine ITERATE. The subroutine 
utilizes the modified Newton-Raphson technique described earlier to solve 
eqn. (5). Note that for combined geometric and material nonlinearity (i.e., 
ANALYS = CNONLIN), the routine GITER is called to obtain a transition non- 
linear solution for the load increment, assuming no additional yielding 
occurs. Then ITERATE proceeds to determine the converged solution which 
includes both geometric and material nonlinearity. The tangential stiffness 
matrix is updated by calling STIFF. For just material nonlinearity (i.e., 
ANALYS = PNONLIN), STIFF is not called. For geometric or combined nonlinear 
analysis, STIFF is called every "NCYCLE” iterations. 

Figure 5 shows a flowchart for the subroutine STRSCAL. STRSCAL calcu- 
lates the incremental stress vector {Aa} corresponding to the calculated 
incremental strains {Ae}. For linear material response, {Aa} is simply the 
product of the constitutive matrix [D] and {be}. For nonlinear material 
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behavior the relationship between {Ac} and {An} depends on the current stress 
state {0^} relative to the yield surface and on the magnitude of the strain 
increment. The relative positions of the stress state and the yield surface 
is determined from eqn. ( 3 ). For convenience in the flowchart, the first term 
in eqn. ( 3 ) is defined to be the effective stress ° e f* F° r an arbitrary 
stress state {a}, the following relationships apply: 

o ({a}) <a * stress state is inside yield surface 
ef ys J 

o r ({o}) =o > stress state is on yield surface 
ef ys 

0 - ( { 0} ) > 0 + stress state is outside yield surface 

ef ys J 

The first step is to calculate the final stress state {o^} assuming no 

additional yielding (block l). Block numbers are indicated at the upper left- 

hand corner of the blocks. If 0 -({0.}) < 0 then {0,} is the correct 

e r I ys 1 

stress state (block 3 A). If not, then {0 } relative to the yield surface is 

o 

examined (block 3 B). If 0 = 0 ({0 }), block 4 B is followed. If 

ys ef o 

a vs ^ a ef^ a o^’ stress state is inside the yield surface. Hence, 

the strain increment must be divided Into two parts: that required to reach 

the yield surface, Ac, and the remainder, Ac, which is the strain increment 

after reaching the yield surface. These strain increments are calculated by 

solving the equations in block 4 A. Next the incremental elasto-plastic matrix 

[D*] is calculated. The final stress state is obtained by adding the linear 

and nonlinear stress increments, [D] {Ac} and [D*] {Ac}, respectively 

(block 6 ). Note that if {0^} had been on the yield surface, {Ac} = 0 and 

{Ac} = Ac. Next the yield stress 0 is updated for strain-hardening 

y s 

materials. Finally, {0^} is scaled back to the new yield surface (block 8 ). 
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Input Data 


The required input data is described in this section* Where applicable, 
the maximum allowable values of the input parameters are noted* 

No. of 

Card set Parameters cards 

1. TITLE(I), I = 1,60 3 

TITLE = TITLE OF PROBLEM 

2. OUTPUT, ANALYS, PLANE, QUADRAT, ENERGY 1 5A8 

OUTPUT = Output option 

= XLONG for long output 

= SHORT for output (the nodal coordinates, element 
connectivity, and boundary conditions are not in the 
output) 

ANALYS = Type of analysis 

= XLINEAR for linear analysis 

= GNONLIN for geometrically nonlinear analysis 

= PNONLIN for materially nonlinear analysis 

= CNONLIN for combined geometric and material nonlinear 
analysis 

PLANE = Option for plane stress/plane strain analysis 
= PSTRESS for plane stress 
= PSTRAIN for plane strain 
QUADRAT - Integration option 

= REDUC for reduced integration 
= XFULL for full integration 

ENERGY = Option for strain-energy release rate calculations 
= DOG for G calculation 
= D0N0 JG for no G calculation 


Format 

20A4 
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Card set 


Parameters 


No. of 
cards 


Format 


3. ITSTEP, NCYCLE, IMAX 1 315 

ITSTEP = Number of steps in the incremental loading 
minimum = 1, maximum = 30 

NCYCLE = Number of iterations between updates of stiffness matrix 
IMAX = Maximum number of iterations allowed before terminating 

4. ACCURACY 1 F10.3 

ACCURACY = Maximum residual allowed in converged solution 

5. NN, NE, NRN 

NN = Number of nodes in the FE model, max. = 650 

NE = Number of elements in the FE model 

NRN = Number of podes with a restrained degree of freedom 

6. Nodal Coordinates: 

x-coordinate 

XX, N(I) = 1,13 * E10.4 , 1315 

XX = coordinate 

N( ) = list of nodes with coordinate XX 

*Input until all x-coordinates are 
specified. End x-coordinate data 
with a blank card. 

y-coordinate 

XX, N(I) , I = 1,13 * E10.4 , 1315 

*Similar to input of x-coordinates 

7. I, IN(I) , JN(I) , KN(I) , LN(I) NE 515 

1,IN,JN,KN,LN = Element number, four node numbers for element I. 

Nodes must be specified in a counterclockwise 
direction. 
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No. of 


Card set Parameters 

cards 

Forma t 

8. K, NRL (2*K-1 ) , NRL (2*K) 

NRN 

315 

K = Node number 




NRL (2*K-1), NRL (2*K) = Constraints in X and Y directions, 

respectively, at node K« 

0 indicates no constraint 

1 indicates constraint 

Note: Do not include degrees of freedom involved in multipoint 

constraints. Do include degrees of freedom with specified 
displacements. 


SKIP 9-12 IF ENERGY = DONOJG 


9. INP 1 15 

INP = Number of node sets used in virtual crack extension 
calculation (maximum = 15) 

10. NEGCAL(I), 1=1, ( INP+1 ) (INP+l)/16 t 1615 

NEGCAL = Element numbers for elements contributing to the nodal 
forces required for virtual crack extension. (See 
example in sketch below. Element numbers are circled. ) 


G> 

(D 

0) 

o 


El 

IE 

MU 

18 

10 

n 

12 

m 

HI 

mm 
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El 
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20^ 

© 

2L~ 

22 


Crack Tip 


■ Crack 


IF INP 

= 

3, 




NEGCAL 

0 

to 

4) 

= 2, 3, 4, 

5 

NFGCAL 

0 

to 

3) 

= 14, 13, 

12 

NDGCAL 

0 

to 

6) 

= 15, 19, 

16, 20, 17, 21 


^Round off to next higher 


integer. 
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Card set 


Parameters 


No. of 
cards 


Format 


11. NFGCAL(I) , I - 1, INF INP/16 t 1615 

NFGCAL(I) = Node numbers for nodes along which virtual 
crack extension forces are calculated. 

List according to distance from crack tip, 
with the crack tip node as the first one. 

(See sketch above.) 

12. NDGCAL(I) , 1=1, (2*INP) 2*INP/16 t 1615 

NDGCAL(I) = Node numbers for the nodes used to calculate 
cracking opening and sliding displacements 


Repeat card sets 13-16 for each material group. 
Maximum number of material groups = 10 
End last group with blank card. 


13. J, XMATER(J) 1 15, A8 

J = Material group number 

XMATER = Material type 

= ELASTIC for linear stress-strain curve 

- ELPLAST for elastic-perfectly plastic stress-strain 
curve 

= BLINEAR for bilinear stress-strain curve 
= RAMOSGO for Ramberg-Osgood stress-strain curve 


14. EX, EY, PYX, GXY 


4E10.3 


E , E , v , G'. 
x> y yx xy 


Young’s modulus in x-direction 

Young’s modulus in y-direction 
£ 

X 

i> = Poisson's ratio 

yx e 

y 

= Contraction in x-direction due to unit 
applied strain in y-direction 


t 


xy 


Shear modulus 


Round off to next highest integer. 
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Card set 


Parameters 


No. of 
cards 


Format 


15. YIELDS, ET, RO, ANM 1 5E10.3 

YIELDS = Yield stress 

ET = Tangent modulus for yielded bilinear material 

RO, ANM = Parameters defining Ramberg-Osgood stress-strain 

ANM 

relation, e = £ + (fj) 

(a) If XMATER = ELASTIC, input YIELDS - ET - RO « 1.0 x 10 21 , 
ANM = 10 

(b) If XMATER = BLINEAR, input proper YIELDS and ET and 
set RO = ANM =0.0 

(c) If XMATER = RAM0SG0, input proper YIELDS, RO and ANM and 
set ET = 0.0 

16. NEL1 , NEL2, NELINC * 315 

NEL1, NEL2, NELINC = Loop parameters used to define 

elements in material group 

NEL1 = First element 

NEL2 = Last element 

NELINC = Loop increment 

e.g. , 1, 50, 20 defines elements 1, 

21, and 31 to be in material group 

*Repeat until all elements in group are defined. 

End card set 16 by specifying NEL1 = NEL2 = NELINC = 0 


17. DELLOAD(I) = 1, ITSTEP ITSTEP/8* 8F10.3 

DELLOAD(I) = Scale factor for proportional load 
vector for load step I. Always 
specify DELLOAD(l) =1.0 


^Round off to next higher integer. 
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Card set 


Parameters 


No. of 
cards 


Format 


18. NLN, NCD, NED 1 315 

NLN = Number of nodes with applied loads 

NCD = Number of multipoint constraints, max = 15 

NED = Number of specified displacements, max = 30 

19. K, FX, FY NLN 15, 2F10.3 

K = Node number 

FX,FY = Loads in x and y directions, respectively 

20. K, KDF, URD NED 215, F10.3 

K = Node number 

KDF = Displacement direction, specify 1 for x direction 

specify 2 for y direction 
URD = Magnitude of displacement 
SKIP 21-24 IF NCD = 0 

21. NMPR(I) , I = 1, NCD NCD/ 16 1615 

NMPR(I) = Number of degrees of freedom involved in the 
Ith multipoint constraint, max = 20 

22. ( (ICDN(I , J) , J=1 ,NMPR(I) ) , 1=1, NCD) NCD sets 1615 

ICDN(I,J) = Jth degree of freedom involved in the 
Ith multipoint constraint 

Start a new card for each multipoint constraint. 

Start with lowest number degree of freedom. 

23. NZKV 1 15 

NZKV = Number of multipoint constraints for 
which there is an applied load 

24. NKV, ATOT NZKV 15, F10.3 

NKV, ATOT: ATOT is the non-zero load associated with 

the NKV set of constrained nodes 
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X 





Figure 2.- Transformed coordinate system for strain-energy-release rate calculation. 


INITIALIZE MATRICES 



STOP 
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READ MATERIAL 
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IF MULTI-POINT 
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ITERATE FOR NONLINEAR 
SOLUTION 

INCREMENT PROPORTIONAL 
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Figure 3.- Flow chart for main program. 
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ITER = 
ITER + 1 



Figure 4.- Flow chart for subroutine Iterate. 
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Figure 5.- Flow chart for subroutine STRSCAL. 
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APPENDIX A 


This appendix gives the names and function of the subroutines and the 
major program variables. 

Subroutines 


NAME 

1. BLKSIGM 

2. BLMAT 

« 3. BMAXQ4 

4. CFILL 

5. CONVERG 

6. DATA 

7. DBAND 

8. DEPMAT 

9. ELPROP 

10. FORCEP 

11. GCAL 

12. GITER 

13. IDVEC 

14. INCLOAD 

15. INITIAL 

16. ITERATE 

17. KLARGE 

18. KSIGNEW 

19. LDATA 

20. LINSOLN 

21. MATMUL 


FUNCTION 

Calculates submatrices for element initial stress matrix 
Calculates nonlinear component of strain-displacement matrix 
Calculates linear component of strain-displacement matrix 
Fills matrix of element nodal coordinates 
Checks for convergence 
Reads nodal coordinate data 

Performs Cholesky decomposition on global stiffness matrix 

Calculates elastic-plastic matrix, [D ] 

ep 

Reads material properties 

Calculates internally generated nodal forces for an element 

Calculates strain-energy release rates 

Solves nonlinear equations 

Fills vector of element degrees of freedom 

Scales load vector 

Initializes variables 

Solves nonlinear equations 

Calculates element large deflection stiffness matrix 

Calculates element initial stress matrix 

Reads load data 

Outputs linear solution 

Performs matrix multiplication 
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22. MULPCON 

23. PROGOPT 

24. RCADD 

25. RESID 

26. RE SUL 

27. SBAND 

28. SDATA 

29. STAXQ4 

30. STIFF 

31. STRSCAL 

32. TRANS 


Modifies stiffness matrix and displacement vector for multi-point 
constraints 

Reads program options 

Adds rows and columns of the stiffness matrix 
Calculates residual force vector 
Calculates strains and stresses 

Solves set of linear equations. (Used with DBAND) 

Reads structural data 

Calculates linear element stiffness matrix 
Assembles global stiffness matrix 

Calculates incremental stresses from incremental strain 
Generates transpose of a matrix 


Program Variables (Arrays are shown with their dimensions.) 


Variable 
AN (1300) 
ANALYS 
ANM 


Definition 

Incremental load vector 
Type of analysis 

Exponent in Ramberg-Osgood equation for uniaxial stress- 
strain curve 


ANTOTAL (1300) 
AR (1300) 

ATOT 

DELLOAD (30) 
DISP (1300) 

DN (1300) 

DPS (30) 

ENERGY 



Total load vector 

Nodal restraint force vector 

Load associated with multipoint constraint 

Scale factor for incremental loads 

Incremental displacement vector 

Total displacement vector 

Vector of specified non-zero displacements 
Option for strain-energy release rate calculation 


25 



FF (1300,4) 

FI (1300,4) 
FXX (10) 

FYY (10) 

IBW 

ICDN (20,15) 
IN (1300) 

JN (1300) 

KN (1300) 

LN (1300) j 
INP 

IPE (1300) 
ITSTEP 
LOADNUM 
MATER (1300) 
MIBW 

MRANK 


Effective stresses at the end of an increment 

Effective stresses at the beginning of an increment 

X-direction forces used in strain-energy-release rate 
calculation 

Y-direction forces used in strain-energy-release rate 
calculation 

Bandwidth of global stiffness matrix 

Degrees of freedom involved in multi-point constraints 

Element connectivity arrays. Connectivity for element 
number I is IN(I), JN(I), KN(I), LN(I) 


Number of node sets used in virtual crack closure calculation 
of strain-energy release rates 

List of yielded elements (only used for output) 

Number of incremental load steps 
Incremental load step number 
Element material group numbers 

Number of columns in global stiffness matrix, SN # 

Currently MIBW =70 

Number of rows in global stiffness matrix, SN. Currently 
MRANK = 1300 


NCD 

NDPS (30) 

NE 

NED 

NLN 

NMPR 

NN 


Number of multipoint constraints 

Vector of degrees of freedom with specified non-zero values 
Number of elements in finite element model 
Number of specified displacements 
Number of nodes with applied forces 

Number of degrees of freedom involved in a set of multipoint 
constraints 

Number of nodes in finite element model 


26 



NND 

Number of degrees of freedom in finite element model before 
applying boundary conditions 

NRL (1300) 

Degree of freedom restraint list 

NRN 

Number of nodes with a restrained degree of freedom 

OUTPUT 

Output option 

PLANE 

Plane stress/plane strain option 

PSI (1300) 

Residual force vector 

QUADRAT 

Integration option 

✓ 

RO 

Parameter in Ramberg-Osgood equation. See definition for 
"ANM" 

SGYBAR (1300,4) 

Current yield stress 

SN (1300,70) 

Global stiffness matrix 

STRESS (1300,12) 

Stresses 

T3 (10,3,3) 

Elasticity matrices for the material groups 

UX (10) 

Tangential displacements vector for strain-energy-release 
rate calculation 

UY (10) 

Opening displacements vector for strain-energy-release rate 
calculation 

X (1300) 

Nodal x-coordinates 

Y (1300) 

Nodal y-coordinates 
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APPENDIX B 


This appendix gives input data and results for three samples problems. 

The first problem (fig. B-la) involves transverse displacement of a long 
thin rod. The finite element mesh is shown with node and element numbers and 
boundary conditions. The left end Is pinned; the right end can move only in 
the "y" direction. The transverse displacement, v, at node 9 was specified 
because the initial transverse stiffness is zero, which would have caused a 
singular stiffness matrix if a transverse load had been specified. Although 
the rod initially has zero transverse stiffness, geometrically nonlinear 
effects stiffen the system as the transverse displacement increases. 

Figure B-lb shows the calculated axial stress in the rod (element 2) as a 
function of lateral displacement. The finite element results are shown as 
symbols. The two curves are exact solutions, derived using simple trigonom- 
etry, for a rod under axial load. One curve is for a linear elastic material 
and the other is for a elasto-perf ectly plastic material with a yield stress 
of 50 KSI. The finite element analysis predicts the nonlinear response very 
well. The differences between the exact results and the finite element 
results are due to the very coarse mesh and the end restraints not being along 
the rod's longitudinal axis. Table B-l lists the numerical values at the 
element centroids calculated by GAMNAS. 

Figure B-2 shows the input data for the linear elastic rod. Required 
changes to this data for the elasto-perfect ly plastic rod are shown in 
parenthesis. 

The second problem involves transverse loading of a double cantilever 
beam. Figure B-3 shows the finite element model, which has 50 nodes and 32 
elements. Two versions of the finite element analysis were used: one version 

used full integration and one used reduced integration. The input data for 


28 



analysis with reduced integration are shown in fig. B-4. The change required 
for full integration is shown in parenthesis. 

The strain energy release rate (using strength of materials) is given by 



A transverse load of 20 lb. was used, resulting in a moment of 40 in. /lb. 

From eqn. (Bl), G is calculated to be 1.92 lb/in. The full and reduced 
integration yielded 1.45 lb/in. and 1.97 lb/in., respectively. Even with a 
coarse mesh, the reduced integration version yielded an accurate result. The 
full integration version illustrates the well-known poor performance of 
isoparametric quadrilaterals in modeling bending deformation. 

The final problem (see fig. B-5a) involves polar symmetric loading of a 
rectangular region. By imposing appropriate boundary conditions along x = 0, 
only half of the region needed to be modeled. The polar symmetric conditions 
are imposed using multi-point constraints to specify u(o,y) = -u(o,-y) and 
v(o,y) = -v(o ,-y) . 

Figure B-5 shows the finite element model before and after loading. 

Table B-2 gives the numerical values of the nodal displacements. The 
required input data are shown in fig. B-6 . 
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TABLE B-l AXIAL STRESS- IN LONG THIN ROD 


LATERAL 

DEFLECTION 

INCHES 

1 

AXIAL STRESS, KSI 
ELEMENT 

2 3 

4 

MATERIAL 

TYPE 

1 

8.868 

8.874 

8.874 

8.868 

LINEAR 

2 

43.48 

43.52 

43.52 

43.48 



3 

104.8 

104.9 

104.9 

104.8 



4 

191.7 

191.9 

191.9 

191.7 



5 

303.7 

303.9 

303.9 

303.7 



1 

8.868 

8.874 

8.874 

8.868 

NONLINEAR 

2 

40.03 

40.09 

40.09 

40.03 



3 

48.61 

48.68 

48.68 

48.61 



4 

49.53 

49.60 

49.60 

49.53 



5 

49.84 

49.91 

49.91 

49.84 




30 













TABLE B-2 NODAL DISPLACEMENTS (INCHES) FOR RECTANGULAR 
REGION WITH POLAR SYMMETRIC LOADS 


NODE 

o 

*— H 

X 

p 

v x 10 4 

NODE 

u x 10 4 

v x 10 4 

i 

-1.095 

.08053 

14 

1.291 

-.4773 

2 

-.6 526 

.02078 

15 

2.146 

-.6216 

3 

.2 x 10 -23 

.2 x 10" 23 

16 

-.4933 

-1.063 

4 

.6 526 

-.02078 

17 

-.4438 

-1.028 

5 

1.095 

-.08053 

18 

.2080 

-1.140 

6 

-.9096 

-.2388 

19 

1.574 

-1.100 

7 

-.4577 

-.2521 

20 

3.415 

-1.268 

8 

.1880 

-.2250 

21 

-.1146 

-.5 x 10" 23 

9 

.9405 

-.2243 

22 

-.3693 

-1.195 

mm 

1.462 

-.3232 

23 

.3345 

-2.047 

n 

-.7632 

-.6573 

24 

1.532 

-2.900 

12 

-.3817 

-.6648 

25 

5.721 

-3.731 

13 

.2632 

-.5427 
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a) FINITE ELEMENT MODEL FOR A LONG THIN ROD 



0 1 2 3 4 5 

Lateral deflection, in. 

b) AXIAL STRESS VS. SPECIFIED TRANSVERSE DISPLACEMENT 


Figure B-l.- Transverse displacement of a long thin rod 








LONt, SLENDER BEAM-COLUMN COSMIC-1 

* # *** ## ****** # * ######################################### 


SHORT 

5 


gnonlin pstress reduc 
20 


DONOJG 


(CHANGE GNONLIN 


TO CNONLIN) 


0 . 010 
10 4 


0 0 1 


5 0 3 

4 

10 . ' 5 

6 

15. 7 

8 

20. 9 

10 

0000000000000000000000000000 

0 . 0 1 

3 5 

■ 2 2 

4 6 

0000000000000000000000000000 

1 1 3 

4 2 

2 3 5 

6 4 

3 5 7 

8 6 

4 7 9 

10 8 

1 1 1 


9 1 1 


1ELASTIC 



0 100E+08 0 1 OOE+OQ 0. 300E+00 0. 385E+07 
0 500E+05 0. OOOE+OO 0. OOOE+OO O OOOE+OO 
14 1 

0 0 0 


(CHANGE ELASTIC TO ELPLAST) 


1 . 000 
0 o 

9 2 


^2. 000 3. 000 4. 000 5 000 

1 . 000 


Figure B 2. Input file for linear elastic rod. Changes required for elasto-plastic 

rod are shown in parentheses. 
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Crack tip detail 



Figure B-3.- Finite-element model for double-cantilever beam. Crack extends from 

X = 0.0 to 2.0 along the line Y = 0.0. 









DOUBLE CANTILEVERED BEAM 


SHORT XLINEAR PSTRESS REDUC 
1 1 1 


COSMIC-3 


DOG 


»**«»**#*#««* ##### 

(CHANGE REDUC TO XFULL) 


1. 000 


50 32 

2 





4000E+00 

1 

2 

3 

4 

5 

2000E+00 

6 

7 

8 

9 

10 

1 OOOE+OO 

1 1 

12 

13 

14 

15 

OOOOE+OO 

16 

17 

18 

19 

20 

1 OOOE+OO 

21 

22 

23 

24 

25 

2000E+00 

27 

28 

29 

30 

31 

4000E+00 

33 

34 

35 

36 

37 

S000E+00 

39 

40 

41 

42 

43 

2000E+01 

45 

46 

47 

48 

49 

1 OOOE+OO 

1 

6 

11 

16 

21 

5000E-01 

2 

7 

12 

17 

22 

OOOOE+OO 

3 

8 

13 

18 

23 

OOOOE+OO 

24 

30 

36 

42 

48 

5000E-01 

4 

9 

14 

19 

25 

1 OOOE+OO 

5 

10 

15 

20 

26 

1 1 

6 

7 

2 



2 6 

11 

12 

7 



3 11 

16 

17 

12 



4 16 

21 

22 

17 



5 21 

27 

28 

22 



6 27 

33 

34 

28 



7 33 

39 

40 

34 



8 39 

45 

46 

40 



9 2 

7 

8 

3 



10 7 

12 

13 

8 



11 12 

17 

18 

13 



12 17 

22 

23 

18 



13 22 

28 

29 

23 



14 28 

34 

35 

29 



15 34 

40 

41 

35 



16 40 

46 

47 

41 



17 3 

8 

9 

4 



18 8 

13 

14 

9 



19 13 

18 

19 

14 



20 18 

24 

25 

19 



21 24 

30 

31 

25 



22 30 

36 

37 

31 



23 36 

42 

43 

37 



24 42 

48 

49 

43 



25 4 

9 

10 

5 



26 9 

14 

15 

10 



27 14 

19 

20 

15 



28 19 

25 

26 

20 



29 25 

31 

32 

26 



30 31 

37 

38 

32 



31 37 

43 

44 

38 



32 43 

49 

50 

44 




20 


26 

32 

38 

44 

50 

27 

28 
29 

31 

32 


33 

39 

45 

34 

40 

46 

35 

41 

47 

37 

43 

49 

38 

44 

50 


45 
50 
1 

19 
18 

23 24 

1 ELASTIC 
0. 100E+08 0. 
0. OOOE+OO 0 
1 32 

0 0 


I00E+08 0. 300E+00 O. 385E+07 
OOOE+OO 0. OOOE+OO 0. OOOE+OO 
1 
0 


1. 000 

1 o o 

50 0. 000 


20. 000 


* Figure B-4.- Input file for reduced integration analysis of double-cantilever beam. 

The change required for full integration is shown in parenthesis. 
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POLAR-SYMMETRIC LOADING OF RECTANGULAR REGION 
***#*****##************************* # ** # ** # ** # 


XLONG XL I NEAR PSTRESS XFULL DONOJG 


1 

1 

1 





1. 

000 






25 

16 

2 





0. OOOOE+OO 

1 

2 

3 

4 

5 

0. 1000E+01 

6 

7 

8 

9 

10 

0. 2000E+01 

11 

12 

13 

14 

15 

0 3000E+01 

16 

17 

18 

19 

20 

0. 4000E+01 

21 

22 

23 

24 

25 

0. OOOOE+OO 

1 

6 

11 

16 

21 

0. 1000E+01 


7 

12 

17 

22 

0. 2000E+01 

3 

8 

13 

18 

23 

0. 3000E+01 

4 

9 

14 

19 

24 

0. 4Q00E+Q1 

5 

10 

15 

20 

25 

1 

1 

6 

7 

2 



2 

6 

11 

12 

7 



3 

11 

16 

17 

12 



4 

16 

21 

22 

17 



5 

2 

7 

8 

3 



6 

7 

12 

13 

8 



7 

12 

17 

18 

13 



B 

17 

22 

23 

18 



9 

3 

8 

9 

4 



10 

8 

13 

14 

9 



11 

13 

18 

19 

14 



12 

18 

23 

24 

19 



13 

4 

9 

10 

5 



14 

9 

14 

15 

10 



15 

14 

19 

20 

15 



16 

19 

24 

25 

20 



3 

1 

1 





21 

0 

1 





1 ELASTIC 







0. 100E+08 0. 100E+08 0. 300E+00 0. 385E+07 
0. OOOE+OO 0. OOOE+OO 0. OOOE+OO 0. OOOE+OO 
1 16 1 

0 0 0 

1. 000 

14 0 

25 1000.000 0.000 
2 2 2 2 

-1 -9 

-3 -7 

-2 -10 

-4 -8 

0 


FILE-COSMIC-5 


Figure B-6.- Input file for polar symmetric loading of rectangular region 



APPENDIX C 


This appendix discusses error messages and potential debug strategies. 

A self-explanatory diagnostic message is output and execution terminated 
under the following conditions: 

1) A node has an unspecified coordinate 

2) An element has an unspecified material group number 

3) The plane stress or plane strain option is spelled incorrectly 

4) An element has a linear stiffness matrix with a diagonal element less 
than or equal to zero. 

5) The rank or bandwidth of the global stiffness matrix exceeds the 
maximum allowed. 

If the global stiffness matrix is singular, the decomposition routine, 
DBAND , prints "Matrix is singular" and halts execution. Failure to specify 
sufficient restraints to prevent rigid body motion is a frequent cause for a 
singular stiffness matrix. A singular stiffness matrix is often encountered 
in geometrically nonlinear analysis because the load increments are too large 
(which causes the iterative solution process to diverge) or because buckling 
occurs. The maximum allowable load increment can only be determined through 
experience. However, frequent updating of the tangential stiffness matrix 
(i.e., a small value is input for NCYCLE) does permit larger load increments. 

The internally generated forces at all nodes are calculated and output. 
These forces should be numerically zero except at nodes where loads are 
applied or displacements are specified, or at nodes involved in a multi-point 
constraint. Errors in modeling will often cause spurious nodal forces, which 
can be used to help isolate the modeling errors. 

Plotting all finite element models is highly recommended, since the plot 
will quickly reveal many input errors. To track down errors not diagnosed by 
GAMNAS , host computer debug utilities are recommended. 
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